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1 Mathematical Definition of Nestedness 



In the main text, we refer to nestedness as the classical NODF measure (nestedness metric 
based on overlap and decreasing fill) [1]. Alternative measures of nestedness are possible. 
In particular, we have also tested our results with the definition of nestedness rj introduced 
by BastoUa and collaborators [2]. In this section we briefiy report and summarize the 
characteristics of these two nestedness measures. Let afj^ be the adjacency matrix of the 
bipartite graph, where the index i (j) runs over the plant (insect) indices: it is 1 if the i-th 
plant is pollinated by the j'-th insect and zero otherwise. The analogous adjacency matrix 
for insects is aff — a^^. Henceforth the set of plant (insect) indices will be denoted 
P (A), whereas the number of plant (insect) species will be denoted by np (n^). The 
number of links (degree) of the i-th plant (insect) is kf = af^^ {kf = afj^) and 
it gives the number of insects (plants) that the i-th plant (insect) interacts with. The 
number of common links (i.e. the number of common pollinators) the i-th and the j-th 
plant have (also called the unweighted overlap matrix) is given by 



oS^^af/af/. (SI) 

k 

In an analogous way, one defines the overlap matrix with respect to insects. 
NODF measure is defined as [1] 



NODF^ 



1 



+ 



A(A-l) 
2 



(S2) 



where T^^ = 0 if kf = k^ , and — '^fjl iiiin(A;f", kf), when both i and j belong to the 
same set X = A, P. 

The alternative definition of nestedness r] is [2] : 

— where r, - (S3) 

All the results presented in this work qualitatively hold for both these measures of nest- 
edness. 

Both nestedness take values in the interval [0, 1], where 1 designates a perfectly nested 
network and 0 indicates a network with no nestedness. They are highly and significantly 
correlated [3]. Prom Eq. (S3) a perfect nestedness is achieved, when the following hierar- 
chy among plants and insects is found: 

Xi C X2 C ... C X„^, for X^P,A (S4) 

where X,- is the subset of species with which species j interacts, and the sub-indices are 
ordered from the most specialist plant (insect), i.e. the species that has the minimum 
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number of inter-species interactions, to the most generalist, i.e. one that has the maximum 
number of inter-species interactions. The main difference between Eq. (S3) and (S2) is 
that, in the latter case, in order to have perfect nestedness, the strict inequahty in Eq. 
(S4) must hold. 

2 Connectance of mutualistic networks 

The connectance Cr is the percentage of mutualistic interactions over all possible ones 
{ua - Up). All the results we have presented in the main text arc obtained for connectance 
that scales with the number of species as Cr ~ A/S^'^, and where S = ua + np. In 
fact, this is the best fit (Table SI) for Ct{S) = aS^ using empirical ecological networks 
(see Figure SI). Data sets of empirical mutualistic networks were taken from the work of 
Rezende et al. [4] .The community size is defined as the total number of plant and animal 
pollinator species. This result is in qualitative agreement with previous results on the 
scaling of the connectivity with the size of mutualistic ecological networks [5] . 

Table SI: Fit results and corresponding goodness of fit tests. 

Cr{S) Estimate Standard Error t-Statistic P- Value 
a 4.21414 1.71557 2.45641 0.0171028 

b 0.776538 0.1107747 7.0101 3.06089-10-9 



For simplicity and in the absence of data, we also set the connectance of the intra- 
species competition interaction matrix in the same way, i.e. Cq ~ 4/5'°'^, but the conclu- 
sions are robust with respect to other choices. 
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Figure SI: Best fit (red solid line) of the connectivity as a function of the number of 
species for 56 mutualistic communities. Dashed gray lines represent the region within 
the ±1 standard deviation confidence interval for the exponent estimate. The plot is in 
log- log scale. 
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3 Mathematical Framework: Community Models 



We present here the theoretical framework and mathematical details of the community 
dynamics model that incorporates mutualistic interactions between plants and pollinators 
or seed dispersers. We start with the general case of a community composed of S species 
comprising np plants and ua animal insects or seed dispersers with S = np + n^- Species 
in the same group are in competition with each other while interacting mutualistically 
with species in the other group. 

x[ and xf denote the abundances of the i-th plant species and the j-th animal polli- 
nator species respectively and x — {xf, xf, ■■■,x^p,x^p-^i, Xg} is a ^'-component vector 
giving the abundance of each of the species in the community. Xi arc variables. The 
intrinsic growth rates of the species in the absence of competition and mutualism is rep- 
resented by the S'-dimensional vector a. Species interact according to a given S x S 
interaction matrix M. The matrix is composed by four blocks, two describing the di- 
rect competitive interactions flpp and Qaa among plants and insects respectively, and 
two assigning the mutualistic interactions TpA and Tap between insects and plants and 
vice-versa. Therefore the interaction matrix M has the following structure: 

d UJi2...UJinp':^inp+l 7l5 \ 

^21 d :72,np+l 

^npl d i7npnp+l ■■■ ■■■ lnp,S (S5) 

7np+l,l 1 d UJnp+lnp+2 ■■■ i^np+lS ' 

7np+2,l ; ■■■ d ...UJnp+2,3 

751,1 :^5np+l d ) 

where Ma = d>0fori = l,...,5' represents the species self -limitation effect. 



Vlpp 


TpA 


Tap 





3.1 Holling Type I Functional Response 

For the sake of mathematical simplicity, we start by representing both direct competi- 
tive and mutualistic interactions between species through a linear functional response of 
HoUing Type I: 

^ = - MijX^ = f{x), (S6) 

The bipartite mutualistic interactions block matrices are built in the following way: for 
each pair 7ij,7jj (with i = l,..,np and j = np + 1, ...S) we draw a random value p 
from a uniform distribution between zero and one (f/[o,i]). If p < Cr, then we set 'jij ~ 
Jji ~ ~|-^(0, CTp)!, otherwise jij = jji = 0. J\f{ii,a'^) is the normal distribution with 
mean /i and variance cr^, and determines the strength of the interactions among species. 
Analogously, the competition interactions ouij and ouji are drawn with probability Cn from 
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|jV(0,(T^)|. We note that because of this construction, the adjacency interaction matrix 
a (where aij = 1 if i is connected to j and zero otherwise), is symmetric. 

The local stability of the underlying population dynamics [6] is governed by the lin- 
earized system of differential equations 6x = ^5x around the stationary point x*, where 
Sx — X — X* and the matrix $ has elements 



_ dfi{x 
(Pij = 



dxj 



-x*M,„ (S7) 



where x* = M~^a and x* > 0 Vi. Therefore the main feature characterizing the local 
stability is the interaction matrix M. 

Let us consider first the case where the stationary point of Eq. (S6) isx* = (1, 1, 1..., 1) 
so that $ = —M and a — M • x*. An approximate, but important indication of the 
stability of the stationary solution x* for the community dynamics given Eq. (S6) can be 
obtained by extending methods recently used to study the stability of structured ecological 
networks [7] to our case (mixture of mutualistic and competition interactions) . It can be 
shown that in the large S limit the maximum real part value of the eigenvalues (A) of the 
linearized system $ is 

<»/x\ J, f {S - 1) {Cr<Jr - Cnaa) Caaa - Crar , . 

maxjR(A) = — a + max I ^ -= , -= h (S8) 

V ^J2^T yl'K 



^ ^fS_{^avC^a^ + (-2Cr + tt + 1) Cpa^ + (-2Cn + tt + 1) Cqct^) 



27r^(2Cr(7rCn(Tn + (tt - Cv) Cvol + (vr - C^) C^al) 

By setting max9?(A) = 0 and solving for (jp we obtain the approximate critical strength 
for mutualistic interactions (for large S) 



^"""^ if(3u<l (S9) 



5(1 - l3u)Cr 



2^d^Cr (tt (/32 + 1,2) -(f3- „)2c^) ^ 
^Cr ((/32 + 1/2) (_2Cr + ^ + 1) + ApuCr) 

+ — ^— — \i^v>\ (SIO 

5((^2 + ^2)(_2Cr + 7r + l) + 4/3i/Cr)' 



where we have set 



Cn = vCt (Sll) 
an = /?ar. (S12) 

Setting 1^ = 1, one has three different situations (see Table S2): if /3 > 1, then competition 
dominates over mutualism, if /3 = 1 competition and mutualism are comparable. Finally, if 
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/9 < 1, the mutualistic interactions dominate the competition. We remark that the results 
presented in Table S2 are only approximate. The extension to the Wigner semicircular 
law ([6]) for an asymmetric real matrix is exact only for the case when the correlation 
between matrix indices is solely between opposite diagonal elements, i.e.{MijMki) = 0 if 
k ^ j and/or I i [8]. For this reason, in order to ensure being within the stability 
region, it is enough to set ar ^ cTc (see Figure S2) 



Table S2: Different scenarios for the interactions strength and corresponding approxi- 
mated stability regions for S* ^ 1. We have imposed that the connectances in the direct 
competition and mutualistic matrices are the same (u = 1). 



Strengths of the interactions 


Scaling Stability Threshold 


weak Competition and strong Mutualism (/3 < 1) 
strong Competition and weak Mutualism (/3 > 1) 
Competition ^ Mutualism (/3 = 1) 


a, ~ 1/{CtS) 
a, ~ l/WCrS) 




Figure S2: Eigenvalues spectrum of random structured ecological network (gray dots) 
with both mutualistic and intra-species competitions interactions as given by Eq. (S7). 
In red are shown the predictions of the extended Wigner semi-circular law [8, 7] of the 
distribution for the eigenvalues of the matrix $ . We have analyzed three possible scenar- 
ios: /3 = 3 where competition is dominant with respect to mutualism, /3 = 0.3 (mutualism 
stronger than intra-species competition), and /3 = 1, i.e. the competition and mutualism 
are comparable. All networks are generated for 5=100, Cr = Cq ~ 0.3, d = 1 and 
CTr = 0.3(Jc. 
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3.2 Rolling Type II Functional Response 

The theoretical framework can be extended to non hnear functional responses. In par- 
ticular, following Thebault and Fontaine [9], we model population community dynamics 
using a HoUing type II functional response: 

^ = (ai - dxi - — ^ = gi{x), (S13) 

at V hij + }_^^dikXkJ 

where dij = 1 ii i E A (P) and j & P {A) and zero otherwise, while hij is a parameter 
having the role of a handling time for the interaction between i and j. In the following, 
for simplicity we set hij — h for all i, j. This choice, although specific, leads to the same 
results one would get with random hij and it allows one to determine the stationary state 
in a much simpler way. We note that in this case we have isolated the self interaction 
term {—dxi) and thus the diagonal elements of M are zero, i.e. Ma — 0 (Vi). Setting 

where xp {xa) is the total plant (animal pollinator) population abundance, the stationary 
solution of Eq. (S13) is given by the self consistent equations 

< = ^[(ci-l + 6(xp,XA))-']ij -a,- (S15) 

3 

Xp = (S16) 

3&P 

XA = J^x*, (S17) 

that can be solved numerically. 

The local stability of the dynamics governed by Eq. (S13) is given by the community 

matrix $ (6x = ^6x) of elements: 

_ dgijxi) _ . Mij + {dx*-ai)eij . 

"^'^ = ^^1^^ - ""^^ ^^'^"^ + h + EjO^,x* J ^^^^^ 

where x* is given by Eqs. (S15)-(S17) and x* > 0 Vi. 
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4 Algorithm for the Maximization of Species Abun- 
dance 



4.1 Community- level Optimization 

In order to implement the optimization principle, the most straightforward dynamics is 
the following. We start with Xi = 1 for i = 1,...,S and a randomly generated matrix 
M as explained above, a is determined as a = Mx* {a = {dl + b)x* if Holling type 
II dynamics is considered). During the optimization dynamics a is kept fixed and the 
matrix entries of F are swapped so to keep fixed the initial entry distribution. At each 
time step, we randomly pick a pair of interacting species and another pair which are 
not interacting (i.e. two pairs {k,l) with 7^ {k,l) and such that 7^ 0 

and 7fc/ = 0), and we switch the interaction between these two elements {i.e.jj^j^^ = 7^^, 
— 7j"). We then calculate the new equilibrium point: x*{n -\- 1) = M^^{n + 1) ■ a. 

If ^j=iX*{n + 1) > Yl,j=i^*j{^)^ then we accept the swap and F = V{n + 1), otherwise 
F = r(n). We repeat this procedure for T steps. We checked that variants of the algorithm 
yields the same results. For instance, if we also allow for non-zero elements to switch, 
the resulting optimal interaction matrix does not change. For a < v^e also performe 
a simulated annealing [10] stochastic algorithm and we did not find relevant differences 
in the final results. Thus the final architecture of M^pt is largely independent of the 
details of the dynamics chosen to perform the maximization of the total population of the 
community. Our results do not change if the starting x* are drawn at random uniformly 
in the range [0.5, 1.5]. 

The total number of time steps T in the optimization procedure is chosen in a manner 
to ensure that in the final steps 90 — 95% of the attempted interaction switches are 
rejected (see Figures S3-S5). The precise number of time steps depends on the particular 
setting for the optimization {S, (Jr, ctq), but for simplicity we have set it to T = 1005". 

4.2 Species-level Optimization 

The above optimization principle may sound as if it is based on some sort of group selection 
mechanism. There is an interesting heated discussion on the role of group/kin selection in 
evolution [11, 12, 13]. Its not our purpose here to enter in this scientific discussion and our 
results do not have to be interpreted as evidence for a particular evolution mechanism. 

In order to not rely on a interpretation of our algorithm based on group selection, we 
performed a variation of the optimization principle: The interactions of the community are 
rearranged in order to increase the abundance of the pollinator or plant species involved in 
the interaction rearrangement. In this way, the switching an interaction improves species' 
fitness (which, in the optimization, is approximated by abundance) and we have a more 
plausible mechanism behind the shift from a random to nested community. 
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0=0.0125 




Figure S3: Number of time steps in order to reach the threshold corresponding to rejection 
of 95% of the attempted interaction switches during the optimization of the community 
population abundance for the case competitionfsmutualism (/3 = 1) and with Cr = Cn = 



The optimization algorithm is similar to the one presented in the previous section: We 
start with Xj = 1 for i = 1, 5* and a randomly generated matrix M; a is also determined 
as before. During the optimization dynamics a is kept fixed and the matrix entries of F 
are swapped so to keep fixed the initial entry distribution. At each time step, we select a 
species, j, randomly and an existing link to one of its partner species k; attempt a rewire 
between the j — k and the j — m links where m is a putative mutualistic partner species, 
that is jjk is interchanged with ■jj^- If the j — m link already exists, i.e. 'jjm is different 
from zero, the switch leads to an interchange of interaction strengths, otherwise the swap 
corresponds to moving the j — k link to j — m. The switch is accepted if and only if it 
does not lead to a decrease of the population abundance of species j in the steady state 
corresponding to the new network configuration. We repeat this procedure for T steps 
(see Figure S6). In order to be sure to reach a local maxima, for this case we set the 
number of steps scale super-linearly with the size T = 100^3/2 

because the population 

dynamics is slower than the previous algorithm. We find that also in this case a nested 
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Figure S4: Number of time steps in order to reach the threshold corresponding to the 
rejection of 95% of the attempted interaction switches during the optimization of the 
community population abundance for the case competition>mutualism (/3 = 3) and with 
Cr = Cn = 4/50-8. 



community architecture emerges for the mutualistic community and all the related results 
still hold. 

4.3 Relation Between Species-level Optimization and Community- 
level Optimization 

The compatibility of individual species fitness and group selection is a convenient fea- 
ture of mutualistic networks. This non-trivial collective behavior can be demonstrated 
in the following way: the stationary community abundance is given by x* = M^-^a*. 
If we impose small variations to the structure of the interactions, i.e. M + 6b, their 
effect propagates to the community populations in a highly non linear way: x* — ?■ 
~ '^jki^ik'^^kiM^^aj. Thus, the local interaction switch has an impact on the 
populations of all the species in the communities, generating a non trivial "collective be- 
havior" (Figure S7), where the benefit to a single species in the long run contributes to 
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a =0.0025 



a =0.0125 




Figure S5: Number of time steps in order to reach the threshold corresponding to the 
rejection of 95% of the attempted interaction switches during the optimization of the 
community population abundance for the case competition<mutualism {(3 = 0.3) and 
with Cr = Cn = 



increase the group fitness (i.e. increase the community population - see Figure S6). 

As we will show in section 5, the relation between species-level optimization and the 
total number of individuals in the community can be also studied analytically using a mean 
field approximation. Indeed we show that an increase in the population of an individual 
species leads to a net increase of the population abundance of the whole community. 

4.4 Optimization while Assembling Communities 

We also tested a more realistic scenario where the interaction network is progressively 
"assembled" over the course of time [14, 15, 16] and meanwhile it is "reorganized" toward 
a more optimal state. We start with four mutualistic communities with five plant and five 
pollinators that randomly interact with probability C. In each one of the communities, we 
maximize the abundance of the individuals in the community as described in section 4.1. 
After 250 attempted swaps we merge communities 1 with 2, and 3 with 4 so to have two 
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Figure S6: Adaptive framework algorithm: at each time step, we randomly pick a species 
i and one of its interacting partners j. We then pick another species k which may or 
may not interact with i. We rewire the interaction between these two elements. If the 
abundance of the initially picked species i does not decrease after the interaction switch, 
then we accept the swap. The parameters used in this simulations are 5* = 20 and 
C = 4/5*°'^. Holling Type 11 dynamics has been used. 



communities with 10 plants and 10 pollinator species. We again turn on the optimization 
for T = 500 time steps. At this point we merge the two communities so to form a 
bigger mutualistic community with S* = 40 species and we continue the optimization for 
T = 2000. As shown by Figure S8, the final architecture is again nested. 

Indeed, our results thus prove to be very robust and largely independent of the details 
of the optimization algorithm. In fact the spontaneous emergence of generalist and spe- 
cialist species from the optimization of the population abundances is independent of the 
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Figure S7: Example of the non trivial collective behavior emergent from the interaction 
switch of a species. The local rewiring between partners i — j and i — I impact the new 
stationary abundances x of all the species in the community. In a single interaction swap, 
some species may be disfavored but in the long run there is an average increase of the 
total population size of both pollinator and plant species (see Figure S6). 



specific dynamical rules of the proposed variational principle as long as stationary condi- 
tions are reached. The dynamical rules are not meant to mimic the actual dynamics of a 
community. Rather, they are efficient means of probing the stationary state characterized 
by the maximal fitness (individual and collective). 



4.5 Swapping algorithm: visualization of the architecture 

In order to provide a visual indicator of the presence of a nested structure in the interaction 
matrix, we have applied the following procedure: 

1. We considered the two off-diagonal blocks of the interaction matrix, containing the 
Plant- Animal Pollinator, TpA {ua x ^p), and Animal Pollinator-Plant, Tap [rip x 
ua), interactions respectively. We applied the next steps to each block separately. 

2. We set all non-zero interactions to one, i.e we consider only unweighted adjacency 

, . PA(AP) PA(AP) .r PA(AP) I n 

matrices 'j^j ' ^ a^j ' = 1 if 7,^- T ^■ 

3. Using TpA we calculate the number of Plant species that interact with each Animal 
Pollinator species, and vice- versa using Tap- For example, for Animal Pollinator 
species i we compute kf = ^"=1 afj^, and for Plant species i we compute k[ = 
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Figure S8: Maximization of the total community population, while assembling commu- 
nities. Parameters used for this simulation are C = 0.35 and cr = <yn = 0.05. If the 
abundance of single populations is optimized, the same results are obtained. 



4. We sort the species according to their kf [kf) values, in decreasing order: (i', j', /c', /') 
with > kp >k§> ... > k^. 

5. We relabel each species following the permutation {i,j,k,...) — )■ {i' , j' , k' , ...), and 
we swap the rows and columns of the original interaction matrix accordingly. 

In order to compare the average structure of many numerical realizations of random 

and optimized interaction matrices with the same S, Cr, riA, np, ar, we applied the 

procedure described above to each of the n realizations, T^''\ where A; = 1, ...,n. Then we 

calculated the average interaction matrix f by averaging each entry over the n realizations: 
f .. _ 1 y^" ^(fc) 

It is important to observe that this method is only one of many possible ways of 
relabeling species, and it may not yield the best visualization of the nestedness of the 
matrix. In order to have a fair comparison, independent of the relabeling procedure 
chosen, what matters is to be consistent and always apply the same procedure. 

5 Analytical Results: Mean Field Approximation 

In the following section, we present a mean field approximation that shows how: a) 
An increase in the population size of one individual species leads to a net increase of 
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the population abundance of the whole community; b) The communities with larger 
total population have an interaction matrix with higher nestedness and vice-versa. The 
approximation consists in considering the 25* x 25* interaction matrix M with constant 
mutualistic/competitive interactions strength: 



M = Mo + y 



1 + n 

0 



0 




" 0 


r" 


1 + 


+ 




0 



where the T superscript denotes the transpose operation and with 



ij 



Cu{l - 5ij) uj > 0 



(S19) 



(S20) 



—7 with probabihty C (7 > 0) 
0 with probability 1 — C 

where C = Cy = Cn is the connectance and Sij the Kronecker S{Sij = 1 if i = j, 0 
otherwise). Therefore, in the mean-field approximation, we have set C ■ S'^ non-zero 
elements of the same magnitude 7 in the mutualistic interaction matrix, while all 
elements are set to Cu in the intra-species competition matrix. 



5.1 Relation between species-level optimization and community- 
level optimization 

Having in mind the population dynamics given by Eq. (S6), the species-level optimization 
algorithm and assuming that a; ~ 7 <IC 1, we now show using perturbative expansions 
that an increase in the population of one individual species leads to a net increase of the 
population abundance of the whole community. We note that the following results still 
hold if 7 and u are random variables drawn from a given distribution (and w ~ 7 ^ 1). 

Consider the interaction rewiring between species k and species m and n as given by 
Figure 9. Then the interaction network changes as M — >■ SM, with 6Mik = 0 (the 
overall total interaction strength is kept fixed). In particular we have: 

{Mkn - Mkm ifi^m 

Mkm-Mkn iii^n (S21) 

0 if i 7^ m, n. 

Similarly, we also have 

{Mnk - Mmk if i = m 

Mmk-Mnk iii^n (S22) 

0 ii i ^ m,n. 

These equations imply that SM^m — SMkn and SM^k — SMnk- Prom Eq. (S6) the 
stationary populations abundances are given by x* — M~^a. After one optimization step 
we thus have 

x* + dx* = {M + 5M)-^a = (1 + M-^5M)-^x*. (S23) 
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Figure S9: Example of swap during the optimization of the population abundance of the 
individual species, presented in section 4.2 



We have assumed that ~ 7 < 1, i.e. M = 1 + e (|| e ||< 1) and thus 6M ~ e. A 
perturbative expansion of Eq. (S23) to lowest order, yields the following change in the 
population abundances 

5x* = -5M-X*. 



From Eqs. (S21)-(S22), we have 

6x1 = - X] ^^kiX^ 



6x1 



6x1 



6x*i 



6Mkm{Xn - 
i 

- 6M^iX* = -6Mmkxl 

i 

0 (/ 7^ n, m, k) 



—6x1 



and therefore 



6x 



tot 



6x1 



(S24) 

(S25) 

(S26) 

(S27) 
(S28) 

(S29) 



We have thus demonstrated that if 6x1 > 0 then 6x1^^ > 0, where fe^^^ is the variation in 
the total number of individuals in the community. The key result is that the optimiza- 
tion of individual fitness leads to an increase of the community fitness (as measured by 
population abundances). 

5.2 Relation between the total abundance of individuals in the 
community and nestedness 

Given Eqs. (S19)-(S21), the inverse of the matrix Mq can be calculated exactly [17]: 



Mo' 



A 0 
0 A 



(S30) 
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with 



1 



l + Cuj{S-2) iii=i 
—Cui if i 7^ J 



a 

-b 



(S31) 



(1 - Cuj)[l + Cuj{S - 1)] 

If we assume that the steady state of Eq. (S6) is x* = (1,1,...,1), then a is fixed to 
Mx*. Without loss of generahty, we can select an initial matrix M such that the vector 
a is constant, i.e. ctj = a, Vi. 

We define the following dynamics for the interaction matrix M: we select two elements 
from r, one interacting (7^- = 7) and one non-interacting (7^; = 0) pair of species, and 
we swap their values, i.e. 7jj = 0 and 7^^ = 7. Using the new matrix M we recalculate 
the populations at the new steady state from the equation x* = M~^d. The total popu- 
lation abundance, x*"*, is then proportional to the sum of all the entries of matrix M~^: 
= '^Sij ^ij^- accept the change if the total population relative to the new M is 
greater than the one calculated with the old M. This dynamics will lead to a matrix M 
that maximizes V- . 

Let us therefore calculate the inverse of an interaction matrix of the form specified in 
Eq. (S19). By observing that M ^ Mq + V ^ (1 + VMq^)Mo, then we have 

M-^ = Mo\l + VMo^)-^ = Mq-i - Mq^VMq^ + Mq-VMo-VMo"^ + ... (S32) 

The sum of all entries of M~^, at the second order approximation (0(7^)), is equal to: 

x'ot = J2m-/ = J2iMo% - J2(^o'VM^% + J2iMo'V M^'VM,-% (S33) 



Inserting Eqs. (S19-S31) in (S33) we have 

s 

E(^o-% = 2j2Aj^2S[a-b{S-l)] 

ij i,j=l 

S 



where we have used 



0 ATA 
AT^A 0 



= 2 J] {AVA)ij = 2[a - b{S - l)]'C,S^(-7) 

ij i,j=l 



(S34) 



kl 



J2iMo'VM,-'VM,-% = E 



ATAT^A 0 
0 AT^ATA 



[a - b{S - 1)]2 J2 [i^Ar^)ki + {r^Ar),i] = (S35) 

kl 

[a - b{S - l)]'[(a + 6)072 - 2b{jCSY], (S36) 
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where 



o — ^ (TkiTkj + TikTjk) 



(S37) 



ij k 



is the overlap of the mutuahstic interaction matrix. 
Thus to second order in 7 we have found 



Jot 



X 



(S38) 
(S39) 

[l-Cuj)[{l + Cuj){S-l)f 

We note that the same conclusions hold for a random competitive interaction matrix Vt 
for which the off-diagonals elements are equal to uj (0) with probability C (1 — C) 



where 
and 



K + Co =^ o (xC ^a;*°* + constant, 
K = 2S[a - b{S - 1)] + 2-fCS^[a - b{S - - b-fCS^) 
C = (a + 6)[a-6(5- 1)]V 



r 



c 




0.10 



Figure SIO: C^^ is the coefficient which relates community abundance Y^- x* to the overlap 
of the mutualistic interaction network (o). For increasing population, we always obtain 
increased overlap, but the correlation decreases as mutualistic (competitive) interactions 
increases (decreases). Here S = 50 and C = Cr = = 4/S'°'^. The value of the 
mutualistic interaction strength is denoted by 7, whereas the competition strength is u 
(mean-field case). 



We have therefore found a direct relationship between the total population of a com- 
munity and the overlap of the matrix representing the species' interactions: 
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Increasing the total population of a community is equivalent to increase the overlap of the 
species' interaction matrix. 

Finally, observing that every definition of nestedness is an increasing function of the plant 
and animal pollinator overlaps, we conclude that: 

Optimizing (increasing) the total abundance of individuals in the community leads to an 
interaction matrix with greater nestedness. 



6 Numerical Results 

6.1 Correlation between the total abundance of individuals in 
the community and nestedness 

The correlation between community abundance and the architecture of mutualistic inter- 
actions revealed by our analytical calculations is also confirmed by extensive numerical 
simulations. Monte Carlo simulations presented in the main text (see Figure 2B) show 
the positive correlation between the total abundance of individuals in the community 
and random matrices characterized by different values of nestedness. In the numerical 
analysis we start with a random interaction M characterized by low NODF. We then 
swap interactions in order to increase the NODF. In this way we have a set of interaction 
matrices characterized by different values of NODF. For each one of them, given a fixed 
a, we calculate the corresponding community population abundance using Eq. (S6) or 
Eq. (S13). We repeat this procedure R times, so to achieve sufficient statistics. Another 
strong piece of evidence for the correlation between community population abundance 
and nestedness is given by the fact that if we minimize the nestedness of the mutual- 
istic interaction network (we accept the swap only if NODF of the mutualistic network 
decreases), then a corresponding decrease in the total population is observed (see Figure 
S12). 

These results demonstrate that population abundance of a mutualistic community and 
the structure of the interactions networks are correlated, independent of the particular 
optimization rule one chooses. 

6.2 Architecture of Optimal MutuaUstic Ecological Networks 

As result of the maximization of the community abundance we obtain a new structure of 
the interaction matrix F, that wc call the "optimal" interaction matrix Topt- In this sec- 
tion, we present the resulting final architecture of Topt for different regimes of interaction 
strengths presented in Table 1. The final optimal shape of the mutuahstic interaction 
matrix F depends on the direct competition matrix Q,. If Q is turned off {ojij = 0 V?,j) 
then, the maximization of the community abundance tends to generate rectangular blocks 
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Figure Sll: Maximization of nestedness of the mutualistic interaction network. We find 
that the corresponding total abundance increases. 
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Figure S12: Minimization of nestedness of the mutuahstic interaction network. We find 
that the corresponding total abundance decreases. 



in the interaction matrix, i.e. there is a subgroup of plants and pollinators in which each 
insect species interacts with all plants and vice- versa (see Figure S14). All nestedness 
results are presented using the NODF measure, but as claimed in section 1 the same 
results hold if the alternative definition {r]) is used, as shown in Figure S13 

When direct competition is turned on (/3 > 0), the optimal mutualistic interaction 
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Figure S13: Comparison between the statistics of nestedness rj - given by Eq. (S3) - 
for random (gray) and optimized (violet) mutualistic networks. We perform: (A) Maxi- 
mization of community population using Rolling Type I dynamics for the case of strong 
competition ctq = 3ar] (B) Optimization of the abundance of the species involved in the 
interaction rewiring for Rolling Type I dynamics and weak competition [an = crr/3); 
(C). Optimization of the abundance of the species involved in the interaction rewiring for 
Rolling Type I dynamics and ctq = ap. All plots are averages over 100 realizations of the 
optimization algorithm with Cq = Cr = 4S'~°'^ and ar < <Jc- 



matrix Topt has an "upper diagonal" like architecture. The particular nested shape of the 
optimal mutualistic interaction matrix depends also on the strength of the mutualistic 
interaction. In fact, as noted in the main text and in SI sec. 5.2, the nestedness of the 
optimized system increases as cxr decreases. Ecologically, this suggests that if mutualistic 
interactions strengths are weak, then it becomes crucial for the community to organize into 
a highly nested architecture in order to maximize its population and thus its persistence. 

Figures S15-S19 show the predicted NODE from our optimization principle for the 
three typical cases of competitive and mutualistic interactions strengths. The simulations 
have been performed using the Rolling Type I saturating functions and for the community- 
level optimization, but similar results are obtained with the Rolling type II saturating 
function or species-level optimization. 

6.3 Null Model Randomizations 

Given that interactions are heterogeneously distributed across nodes, the question then 
arises whether they are distributed randomly or not. Note that it is not meaningful to 
compare absolute values of nestedness across networks of different sizes, connectance, or 
with different degree distributions[18, 19]. In order to address this issue, we will follow 
previous works [5, 9, 18] and define two null models to benchmark the degree of nestedness 
(see Figure S15). 

In Null Model 0 (also indicated as "random"), we keep fixed the size and the connec- 
tivity of the optimized networks, while the edges are placed randomly within the matrix 
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Figure S14: Architecture of a random mutualistic network Tran and its evolution driven 
by the maximization of population abundance principle. F at intermediate time is in- 
dependent of the presence or absence of direct competition Q, while the final optimal 
network structure does depend on Q. The second and the third interaction matrices from 
the left are averages performed over several runs of the iteration procedure and the scale 
of grey is a measure of the probability that a given entry is present. Time T is measured 
in terms of the number of iterations. All networks are generated for S=50, C ~ S^'^ and 
a <^ (Tc. 



[5]. Moreover, since nestedness changes along with changes in the degree distribution 
[18, 19], we also perform a more stringent randomization (Null Model 1), where we 
keep fixed the average number of connections for each plant and insect, i.e., each element 
A[j of the bipartite randomized adjacency matrix is one with probability kf ■ kf/L and 
0 otherwise {kf, k^ are the degree of plant i and insect j, respectively and L is the total 
number of mutualistic links in the network) [5]. For each optimized network we compute 
100 realizations of null model 1 and then calculate the corresponding average nestedness 
NODFnmi- In this way we can also define the relative nestedness NODF* [9], i.e., the 
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Figure S15: Histograms of the nestedness probability density (PDF) for randomizations 
under null model 0 (gray), under null model 1 (in green) and mutualistic networks, where 
community population has been maximized using Rolling Type I (in red) and Type II (in 
violet) population dynamics. All plots are averages over 100 realizations of the abundance 
maximization algorithm for = and Cq = Cr = 45'~°''^. 



nestedness normalized to null model 1 (nml): 

NODF' = ^ODF -NODF„„,, 
iVOCF„„„ 

As we can see from Figure S16 and Figure 3 in the main text, the relative nestedness 
increases along with the optimization of the interaction network. 
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Figure S16: Relative nestedness as a function of the number of optimization steps T for 
the case of community-level optimization and using Holhng Type II saturating function. 
We note that the results depend on the optimization dynamics. Parameters used in this 
simulation are S=50, Ccamma = 4S'~°'^, Cq = 0 and ar = 0.15. 
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Figure SI 7: Nested architecture measured by NODF versus the number of species of 
the optimized mutuahstic interaction network for three different values of interaction 
strengths for the case competition= mutuahsm (/3 = 1). Community population has been 
maximized using HTI saturating function. The dashed lines give the 1 standard deviation 
confidence interval. 
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Figure S18: Nested architecture measured by NODF versus the number of species of 
the optimized mutuahstic interaction network for three different values of interaction 
strengths for the case competition> mutuahsm (/3 = 3). Community population has been 
maximized using HTI saturating function. The standard deviation confidence interval 
behaves qualitatlvly as in Figure SI 7. 
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Figure S19: Nested architecture measured by NODF versus the number of species of 
the optimized mutuahstic interaction network for three different values of interaction 
strengths for the case competition< mutuahsm (/3 = 0.3). Community population has 
been maximized using HTI saturating function. The standard deviation confidence interval 
behaves qualitatlvly as in Figure SI 7. 
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7 Stability Analysis of the Optimized Networks. 



Community persistence and stability are known as important dynamical properties char- 
acterizing ecological networks, but the way in which the two are related in real systems is 
far from trivial. Community persistence has been considered as a possible driving force 
sculpting network architectures, i.e. the observed life web structures are the ones leading 
to high community persistence. Persistence can be defined as the fraction of initial species 
remaining once the dynamics has driven the system to a stable state. Consequently it 
is relevant only for unstable systems, and moreover, it is sensitive to initial conditions 
and transients, and to the distance of the system from stationarity. It is hard to disen- 
tangle the role of distinct ecological variables on the networks topology in a dynamically 
evolving system, since they all typically vary simultaneously. For instance, the reduction 
of the effective ecosystem diversity due to species extinction that occurs while reaching 
stationarity, necessarily increases the nestedness of the corresponding interaction matrix 
because, on average, even randomly assigned matrices of smaller size are more nested [20]. 
Therefore here we focus on the study of community stability for the optimized stationary 
networks. 

7.1 Analytical Results: relation between Rarest Species and 
Community Resilience 

Using perturbative expansion techniques, we find that rare species control community 
stability. In particular, we find that the dominant real part of the eigenvalues of $ is 
related to the population of the rarest species through 



when (Tr <^ 1. 

7.1.1 Holling Type I 

From Eq. (S7), we know that the stationary populations of the species x* > 0 are 
related to the linearized matrix $ of the Holling Type I dynamics Xi = fi{x) through 
0ij = —x*Mij. The interaction matrix M is built so that 



where the index X = F, f2 depending on whether we are considering the mutualistic or 
competitive block and S is the Dirac 5— function. Thus Eq. (S7) can be rewritten as 



Max[Re{X)] 



Min[x\, 



(S42) 




(f)ij = -x*Mij = -x*Sij - (1 - Sij)Mijx* ^(j>l + V, 



(S43) 
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where 0?- = —x*Sij and Vij — —{1 — 6ij)MijX*. The eigenvalues of 0° are —x*, while its 
eigenvectors are u^'^\0, 0, 0, 1, 0, ...0) where the 1 is the i—th component of u^'^\ We can 
thus calculate the eigenvalues of $ using perturbation theory. If Aj is the i—th eigenvalue 
of $, then we have (cr = maXx(Jx) 

from which Eq. (S42) is achieved. We have numerically tested the prediction given by 
Eq. (S44) and the result is shown in Fig. S20. 

7.1.2 Holling Type II 

This result is also valid for a community dynamics with Holling Type II saturating func- 
tion. In fact, from Eq. (S18) we have that 

h + E,- OijX] 

where this time = -dx*6ij and Vij = -x* ^^^hj^^'g.^J^'^ ■ 

The eigenvalues of are —dx*, while its eigenvectors are u^'^^O, 0, 0, 1, 0, ...0) where 
the 1 is the i—th component of u^'^\ Therefore also in this case we have that the first term 
in the perturbative expansion goes to zero, i.e. u^^^ ■Vu^'^'* = 0 (remember that Ma — 0 and 
9ii — 0 Vi). Consequently, given Aj as the i—th eigenvalue of we have (cr = maXx(Tx) 

Xi = -dx* + o(a^) (S46) 



5jj _ 



+ \}'^'a.y ' ] =^3 + yr,. (S45) 



7.2 Numerical Results: Spectrum of the Optimized Networks 

We performed an analysis of the eigenvalues A of the matrix $ associated with the optimal 
interaction matrix M^t corresponding to the maximal population Xopt (see Eq. (S7)) in 
order to study the stability of the stationary states of the optimized networks. 

Figures S21-S23 show the value of the maximum real part of the eigenvalues of $ after 
the optimization procedure for the three standard cases when Eq. (S6) is used for the 
population dynamics. We find that, after the optimization, the system remains stable 
with high probability if the starting random interaction matrix is stable (dr < a"c), while 
if the initial system is critical, i.e. at the edge of stability (err ~ cr^) or unstable {a > ac), 
then also the optimized community is critical or unstable, respectively. At the same time 
we observe that there is a shift of the optimized system towards the edge of stability, i.e. 
Max(i?e(A)) gets closer to zero. 
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Figure S20: Numerical analysis of the relationship between the community resilience 
(given by the maximum real part eigenvalue of the interaction matrix $) and the pop- 
ulation of the rarest species for different values of ar and for the case competition= 
mutualism (/3 = 1). Red squares represent the real part of the dominant eigenvalue of 
$ for different realizations of the optimization starting from a random interaction ma- 
trix. The linear correlation (black line) predicted by our analytical calculation based on 
perturbative techniques holds within the stability region ar < <Jc < 0.2. 
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Figure S21: Maximum real part of the dominant eigenvalue of $ before (blue squares) 
and after (red circles) the optimization procedure for the case competition> mutualism 
(/3 = 3). 
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Figure S22: Maximum real part of the dominant eigenvalue of $ before (blue squares) 
and after (red circles) the optimization procedure for the case competition< mutualism 
(/3 = 0.3). 
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Figure S23: Maximum real part of the dominant eigenvalue of $ before (blue squares) 
and after (red circles) the optimization procedure for the case competition?^ mutualism 

(/3 = !)• 
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Figure S24 shows how also for Holhng Type II dynamics, the optimized mutualistic 
community is less stable than the random starting one. In fact, another way to check 
the stability of the optimized community matrix is to study the population of the rarest 
species in the mutualistic community, as shown in previous section. 




0 10 20 30 40 50 60 70 80 90 100 
Realizations 

Figure S24: Population size of the rarest species before (dashed gray line) and after (red 
solid line) the single species optimization using Holling Type II saturating function for 
the population dynamics. The x axis represent different realizations of the optimization 
starting from the same initial condition with = (1,1,1...,!). The blue solid line likewise 
shows the population of the most abundant species in the optimized community in differ- 
ent realizations. This spontaneous symmetry breaking between species is a hallmark of 
nested architecture. Parameters of the simulations are S = 20, C = 4:S~^'^, ar = ctq = 0.1 
and h = 0.1. 
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